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Aggregation of rod-like polyelectrolyte chains in the presence of monovalent counterions 
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The Institute of Mathematical Sciences, C.I.T. Campus, Taramani, Chennai 600113, 
India 

(Dated: 3 March 2013) 

Using molecular dynamics simulations, it is demonstrated that monovalent counterions can induce aggrega- 
tion of similarly charged rod-like polyelectrolyte chains. The critical value of the linear charge density for 
aggregation is shown to be close to the critical value for the extended-collapsed transition of a single flex- 
ible polyelectrolyte chain, and decreases with increasing valency of the counterions. The potential of mean 
force along the center of mass reaction coordinate between two similarly charged rod-like polyelectrolytes is 
shown to develop an attractive well for large linear charge densities. In the attractive regime, the the angular 
distribution of the condensed counterions is no longer isotropic. 



I. INTRODUCTION 

Polyelectrolytes (PEs) are charged polymers in a so- 
lution containing neutralizing counterions 1-3 . They are 
common in biological systems, examples include DNA 4 , 
F-actin and microtubules^. A special case is a rod-like 
PE (RLPE) whose persistence length is of the order of 
its contour length. RLPEs are of interest theoretically 
because their thermodynamics is similar to the theoreti- 
cally well studied idealized system of a charged cylinder 
with neutralizing counterions. 

PEs, even if similarly charged, may attract each other 
in the presence of counterions. The aggregation of 
similarly charged PEs has been extensively studied ex- 
perimentally^—, theoretically-^— and numerically^ 6 - - — . 
Despite these studies, the role of counterion valency in 
aggregation remains unclear. Though it has been shown 
unambiguously that multivalent counterions can mediate 
aggregation^^— i 26-^o i 32-^5 j it is still being debated 
whether monovalent counterions can cause aggregation 
in the absence of multivalent salts. There are some ex- 
perimental^— and theoretica l 18 ' 24 results that argue for 
monovalent counterion induced aggregation of PEs. At 
the same time, other experimenta l 11 ' 13 , theoretica l 17 ' 36 ' 37 
and numerical 26-31 studies argue or report the absence of 
aggregation in the presence of monovalent counterions. 

There are different proposals for understanding the 
attraction between two similarly charged PEs. One of 
them suggests that the attraction is induced by the cor- 
related longitudinal fluctuations of the condensed coun- 
terion density^ 6 -. This theory predicts attraction in the 
presence of multivalent counterions, in agreement with 
early numerical simulations 26 , for a range of system pa- 
rameters. The validity of the Gaussian approximation 
made in the theory was questione d 36 ' 38 , and more realis- 
tic, but simple models were considere d 17 ' 36 . These mod- 
els assume localization of the condensed counterions on a 
finite number of allowed sites around the PEs, and were 
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able to explain multivalent counterion induced attraction 
in terms of spatial distribution of the condensed counte- 
rions around the PEs. However, these theories predict a 
very weabi 7 - or zero^ 8 - attraction for the case of monova- 
lent counterions. 

An alternate approach that explains the attraction of 
PEs is by Manning et. al-A 8 -, and is based on the classical 
condensation theory 3 ^. Though this approach does not 
explicitly take into account correlations between coun- 
terions, it predicts attraction between similarly charged 
PEs. In fact, the theory predicts stronger attraction in 
the presence of monovalent counterions than divalent or 
trivalent counterions. The origin of attraction is similar 
to that of covalent bonds, in that the condensed coun- 
terions are shared by the PEs. This theory was later 
extende d 23 ' 24 , using extended condensation theory 4 ^, to 
calculate the interaction free energy as a function of the 
separation between the PEs, linear charge density of the 
PEs and valency of the counterions. The analysis of the 
free energy shows that the attraction is a consequence 
of the increase in the number of condensed counterions 
and the counterion condensation volume with decrease 
in the separation between the PEs. This approach dif- 
fers from the Manning theory in that the assumption of 
infinite dilution of added salt is absent. Although these 
theories predict the formation of stable aggregates, medi- 
ated by monovalent counterions, a clear experimental or 
computational confirmation is still lacking. Most earlier 
numerical simulations 2 ^ 31 ' 34 do not observe attraction 
that is strong enough to form stable aggregates. 

A different phenomena in PE systems, mediated by 
counterions, is the extended-collapsed transition of a 
single flexible^— or semi-flexible PE 47 . As the lin- 
ear charge density exceeds a critical value, an extended 
flexible PE chain undergoes a first-order transitio n 42 ' 44 
to a globular phase, while an extended semi-flexible PE 
transforms into a toroid-like phase 47 . As for aggregation 
of PEs, there are different proposals for the underlying 
mechanism of the extended-collapsed transition. In one 
of the approaches, a one-component plasma model of the 
PE system was shown to have negative compressibility 
when the linear charge density of the PE chain exceeds 
a critical valued. The negative compressibility causes 
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the PE chain to collapse, with the transition being first- 
order. In an alternate approach^ 3 -, the transition was 
studied using a two-state model, where the free energy 
of the extended and collapsed phases of PE chain was 
calculated. The collapsed phase was approximated by an 
amorphous ionic solid made up of the condensed counte- 
rions and the monomers of the PE chain, and the strong 
correlation between the counterions and the monomers 
was argued to be responsible for the collapse of the chain. 
In an another numerical study 41 -, the condensed counteri- 
ons were shown to form dipoles with the monomers of the 
chain, and the resultant dipole-dipole attraction causes 
the chain to collapse. A different approach^, based on 
the dependence of the free energy of a semi-flexible PE 
chain on counterion density fluctuations, argued that the 
extended phase of a semi-flexible PE chain destabilizes 
into a toroid-like phase at high linear charge densities. 

In this paper, using molecular dynamics simulations 
(details in Sec. [TTJ) , we demonstrate the aggregation of 
similarly charged RLPEs in the presence of monovalent 
counterions when the linear charge density of the polymer 
backbone is higher than a critical value (Sec. IIII A|) . We 
argue that the critical backbone charge densities required 
for the onset of aggregation of RLPEs and the extended- 
collapsed transition of a single flexible PE chain are 
nearly equal (Sec. IIII Bj) . We also measure the poten- 
tial of mean force between two RLPEs along the dis- 
tance between their centers of mass and show the exis- 
tence of an attractive well for large linear charge densi- 
ties fSec. IHICj) . In this attractive regime, a spatial rear- 
rangement of counterions around a RLPE occurs and is 
quantified in Sec. IIII D] Sec. IIVI contains a summary and 
discussion. 



II. MODEL AND SIMULATION METHOD 

We model RLPE and flexible PE chains as N spheres 
(monomers), each with charge +qe, connected through 
springs. The counterions are modelled as spheres with 
charge —Zqe, where Z = 1,2 and 3 for monovalent, di- 
valent and trivalent counterions respectively. All counte- 
rions have same valency, and the number of counterions 
is such that the system is overall charge neutral. The 
polymer chain and the counterions are assumed to be in 
a medium of uniform dielectric constant e. 

The interactions between particles i and j are of four 
types: 

Coulomb interaction: The electrostatic energy is given 

by 



U c {ri 



(1) 



where nj is the distance between particle i and j, and qi 
is the charge of the i th particle, and eo is the permittivity 
of free space. 

Excluded volume interaction: The excluded volume in- 
teractions are modelled by the Lennard- Jones potential, 



ULj(nj) = 4ej 
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where is the minimum of the potential and a is the 
inter-particle distance at which the potential becomes 
zero. We use reduced units, in which the energy and 
length scales are specified in units of ey and a respec- 
tively. The depth of the attractive potential and 
its range a are set to 1.0 for all pairs of particles. We 
use shifted Lennard- Jones potential in which Ujjjfrij) 
smoothly goes to zero beyond a cut off distance r c . The 
value of r c is chosen to be 1.0 such that the excluded 
volume interaction is purely repulsive for all pairs, mim- 
icking polymers in good solvents. 

Bond stretching interaction: The bond stretching en- 
ergy for pairs that are connected through springs is given 
by 



(3) 



where kb is the spring constant and b is the equilibrium 
bond length. The values of kb and b are taken as 500 and 
1.12 respectively. This value of b is close to the minimum 
of Lennard-Jones potential, and the spring constant kb is 
large enough so that the bond length does not change 
appreciably from b. 

Bond bending interaction: The rigidity of the polymer 
backbone is controlled by a three-body interaction given 
by 



U e (6) = k g [l + cos( 



(4) 



where 9 is the angle between two adjacent bonds. For 
RLPEs, we choose kg = 10 3 to ensure rigidity of the 
backbone, while kg is set to zero for flexible PEs. 

The relative strength of the electrostatic interaction is 
parameterized by a dimensionless quantity A: 



A 



q 2 £ B 



where £b is the Bjerrum length , 
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where kg is the Boltzmann constant and T is tempera- 
ture. In our simulations, we vary A from 0.22 to 10.93. 

The equations of motion are integrated in time 
using the molecular dynamics simulation package 
LAMMPS 49,50 . The simulations are carried out at con- 
stant temperature (T=1.0), maintained through a Nose- 
Hoover thermostat (coupling constant = 0.1) 51 i 52 . The 
system is placed in a cubic box with periodic boundary 
conditions. We use the particle-particle/particle- mesh 
(PPPM) technique 53 to evaluate the energy and forces 
due to the long range Coulomb interactions. The time 
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step for the integration of the equations of motion is cho- 
sen as 0.001. 

In the current study, three systems are considered. 
In the first system, aggregation of 50 similarly charged 
RLPEs with varying linear charge densities neutralized 
by counterions of varying valency is studied. In the sec- 
ond set of simulations, we locate the extended-collapsed 
transition point of a single flexible PE chain in the pres- 
ence of counterions of different valency. In the third set 
of simulations, we measure the potential of mean force 
between two RLPEs in the presence of monovalent coun- 
terions to understand the observed aggregation. 



III. RESULTS AND DISCUSSION 

A. Aggregation of rod-like polyelectrolyte 
chains 

To study aggregation, we consider a collection of 50 
RLPE chains of 30 monomers each. The system is charge 
neutralized with monovalent, divalent or trivalent counte- 
rions. The density of the system is chosen to be 4.4 x 10~ 4 
monomers /a 3 . At this density, the mean separation be- 
tween the chains is much larger than the length of the 
chains. We vary the non-dimensional parameter A from 
0.22 to 10.93. 

We first show that monovalent counterions can medi- 
ate aggregation of similarly charged RLPEs, in contrast 
to the previous experimenta l 11 ' 13 , theoretica l 17 ' 36 ' 37 and 
numerical^ 6 - - — studies. In Fig. QJA)— (D), we show snap 
shots of the system at increasing times for A = 9.43, 
when all counterions are monovalent. For this value of 
A, electrostatic interactions are much larger than ther- 
mal energies, and aggregation of RLPEs can clearly be 
seen. 

We quantify aggregation by calculating the dependence 
of average aggregate size on A. Two PE chains are de- 
fined to form an aggregate if the distance between any 
monomer of the first chain and any monomer of the sec- 
ond chain is less than 2a. Similarly, a PE chain is defined 
to be part of an aggregate of size m (m > 2) if the dis- 
tance between any of its monomers and any monomer of 
the other m — 1 chains is less than 2a. Other similar defi- 
nitions for what constitutes an aggregate can be found in 
the literatures 7 .! 3 ^, but we find that the results do not 
change noticeably with different choices of the definition. 

Let N m be the number of aggregates of size 
m. Then, the average aggregate size is given by 



Era=i m *™/E m =i' l 'ra' Measuring the equilibrium 
value of N m turns out to be difficult in the aggregate 
regime because of large equilibration times and limited 
computational time. We, therefore, measure N m after 
starting from two different types of initial conditions and 
discarding a fixed number of initial simulation steps. In 
the first type, we start from an initial random configura- 
tion of RLPEs and discard the first 3 x 10 7 time steps. 
We then average the aggregate size over the next 10 7 
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FIG. 1. The two dimensional projection of the RLPEs at time 
steps (A) 0, (B) 1.90 x 10 6 , (C) 8.90 x 10 6 and (D) 4.72 x 
10 7 . The counterions are monovalent and are not shown for 
clarity. The RLPEs aggregate, with the number of aggregates 
decreasing in time. The data are for A — 9.43. 



time steps. The average aggregate size, thus computed, 
is shown by crosses in Fig. [2JA) as a function of A for 
counterions of different valency. The aggregation begins 
at critical values A — 6.57, 2.50 and 1.45 for monovalent, 
divalent and trivalent counterions respectively, which cor- 
respond to an average aggregate size of 2.0. At high 
enough values of A, all the 50 chains form a single aggre- 
gate for both divalent and trivalent counterions. How- 
ever, in the case of monovalent counterions, even for the 
highest value of A (10.93) that we studied, we did not 
observe formation of a single aggregate within the maxi- 
mum duration of our simulations (5 x 10 ), and the final 
configuration of the system consists three aggregates of 
approximately 16 PEs each. In the second set of simu- 
lations, we start with a high value of A and an initial 
condition where all the chains are in an aggregate. We 
then discard the first 10 7 time steps, and measure the 
aggregate size over the next 10 6 time steps. A is then 
decreased and the initial condition is chosen to be final 
state of the previous A value. The data thus obtained are 
shown by circles in Fig.^A). The values for A for which 
the two kinds of runs give the same value for average ag- 
gregate size can be taken to be the correct equilibrium 
value. The formation of a hysteresis loop is clearly seen. 

The difficulty in equilibration in the aggregating 
regime is due to large diffusion time scales. As small 
aggregates form, they take longer and longer to diffuse 
and come close to each other before the next aggregation 
event can occur. We find that for values of A larger than 
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FIG. 2. (A) The variation of the average aggregate size with 
A. The data shown by crosses are obtained from simulations 
where the RLPEs are initially randomly distributed. For the 
data shown by circles we start at high value of A where ini- 
tially all RLPEs are part of a single aggregate. The initial 
configuration for a lower value of A is the final configuration 
of the previous value of A. (B) The relative fluctuation Xc 
in electrostatic energy of a single flexible PE as defined in 
Eq. (J7J. The height of the peaks have been rescaled for clar- 
ity. (C) The variation of asphericity Y [see Eq. (JSJ) ] with A of 
a single flexible PE chain. 



the critical values, the aggregation kinetics is indepen- 
dent of A and the valency of the counterions. In Fig. [31 
we show the time evolution of the history averaged num- 
ber of aggregates (N a ). At the start of the simulations, 
the PEs are randomly distributed ((N a ) = 50). After 
initial transients, corresponding to counterion condensa- 
tion, (N a ) decays as a power law t~ T , with r w 0.66, 
independent of valency. If we model aggregation as irre- 
versible coalescence of point-sized particles (see Refj 54 i 55 
for a review), then the above power law decay corre- 
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FIG. 3. The variation of {N a ), the average number of ag- 
gregates, with time t for different values of A and valency. 
The data is log-binned and averaged over five different initial 
configurations. The solid black line is a power law ~ i _2//3 



sponds to the diffusion constant of an aggregate of size 
m being proportional to m -1 / 2 . 

In the earlier simulations 2 ^—, the interaction between 
similarly charged RLPEs in the presence of monovalent 
counterions was found to be repulsive. Typical values 
of A used in these simulations are 2.1Q2&, 2.90 2 -L 2 2 and 
4.17 28 , which are much smaller than the critical value 
(« 6.57) above which we observe the aggregate forma- 
tion. We suggest that the absence of attractive interac- 
tions required for aggregation in these simulations can be 
attributed to the inappropriate values of A. 



B. Extended-Collapsed transition of a single 
flexible polyelectrolyte chain 

In our earlier simulation s 44 ^ 45 , we had demonstrated 
that a single flexible PE chain undergoes a first-order 
transition from an extended to a collapsed configuration, 
mediated by the counterions, when the parameter A ex- 
ceeds a certain critical value. To establish a possible 
relation between the aggregation phenomena of RLPE 
chains and the extended-collapsed transition of a sin- 
gle flexible PE chain, we study the extended-collapsed 
transition of a single flexible PE chain of 600 monomers 
in the presence of monovalent and multivalent counteri- 
ons. The density of the system is chosen as 2.7 x 10~ 6 
monomers/cr 3 such that the direct contact between the 
PE chain and it's periodic images is avoided. The initial 
configuration of the PE chain is randomly chosen and the 
counterions are uniformly distributed inside the simula- 
tion box. The system is allowed to equilibrate for 10 7 
steps and the averages are taken over a production run 
of 10 7 steps. 

A useful quantity to study the extended-collapsed 
transition is the electrostatic energy per monomer E c . 
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E c shows different scaling with the number of monomers 
in the extended and collapsed phases^. The relative fluc- 
tuation Xcin E c is defined as 



Xc 



N [{El) - {E c f 



(7) 



Fig. HJB) shows the variation of Xc with A. Xc has a 
peak around 7.5,2.57 and 1.40 for monovalent, divalent 
and trivalent counterions respectively. These peaks cor- 
respond to the extended-collapsed transition of the single 
PE chain. 

The extended-collapsed transition can be further 
quantified by studying the variation of asphericity of the 
PE chain as a function of A. We define asphericity as 



Y = 



Ai 



A2+A3 
2 



Ai + A2 + A3 



(8) 



where Ai^.3 are the eigenvalues of the moment of inertia 
tensor with Ai being the largest eigenvalue. The moment 
of inertia tensor G is 



1 N 

i=l 



(9) 



where ri a is the a th component of the position vector 
fi . Asphericity Y is zero for a sphere (collapsed globule) 
and one for a linear rod (extended configuration). For 
all other configurations, it has a value between zero and 
one. 

The variation of Y with A is shown in Fig. [2jC). For 
all the three valencies, Y makes a transition from a value 
close to one (extended configuration) to zero (collapsed 
configuration) at a critical value of A. The transition in 
Y occurs at the same critical value of A where Xc has a 
peak, and corresponds to the extended-collapsed transi- 
tion of the PE chain. We note that these critical values 
roughly coincide with the critical values for the aggrega- 
tion of RLPE chains in the presence counterions of the 
corresponding valency. This suggests the possibility that 
the underlying mechanisms of the collapse of a single PE 
chain and the aggregation of RLPE chains are closely 
related. 




FIG. 4. The reversible work W needed to bring two rod-like 
PEs from a separation of 5<r to d, in the presence of mono- 
valent counterions. The averages are taken over a production 
run of 10 7 steps, after the system has been equilibrated over 
10 7 steps. 



short-ranged, and hence W(d) can be approximated by 
the work done in bringing the rods to a separation d 
from a finite separation d', where d' > d. W(d) is then 
J^dxf(x), where f(x) is the normal component of the 
force required to keep the two rods at a separation x. 



For the evaluation of W(d), we consider two parallel 
RLPEs of 30 monomers each. The system is charge neu- 
tralized by monovalent counterions. The density of the 
system is chosen to be 4.4 x 10~ 4 monomers/u 3 . Fig. 0] 
shows the variation of W(d) per monomer for different 
values of A. For low values of A, W(d) is always pos- 
itive and decreasing with d, showing repulsion between 
the rods. At high values of A, W(d) develops a mini- 
mum which is the onset of attraction between the rods. 
When the depth of the potential s comparable to thermal 
energy T, the attraction will be strong enough to form 
stable aggregates. 



C. The effective interaction potential between 
two rod-like PE chains with monovalent 
counterions 

To understand the observed aggregation of RLPEs 
in the presence of monovalent counterions, we measure 
the effective interaction potential between two similarly 
charged RLPEs. The effective interaction potential, for a 
separation d between the rods, is equal to the reversible 
work W(d) done in bringing them from infinite separa- 
tion to a separation of d. Since the system is overall 
charge neutral, we expect the interaction potential to be 



A similar evaluation of W(d) for two similarly PEs in 
the presence of monovalent counterions was carried out 
in earlier simulation a 29 ' 30 . These simulations were per- 
formed at a range of values of A which is much less than 
the values at which we observe the attractive part in 
W(d), and these simulations failed to observe the attrac- 
tive part in W(d). In RefJ^, solving a simple model of 
PE system with monovalent counterions, the free energy 
was evaluated as a function of the separation between the 
rods. For the value of A studied (A = 4.10), it was shown 
that the depth of the attractive well was insufficient to 
bind the two PEs. 
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D. The angular distribution of the condensed 
counterions 

It has been observed in experiments^ 6 - and simula- 
tions^ that the angular distribution of the condensed 
multivalent counterions changes as the PEs approach 
each other. Many theoretical studie d 17 ' 36 ! 38 also argue 
that the spatial arrangement of the condensed multiva- 
lent counterions around the PEs plays an important role 
in developing an attractive interaction between the PEs. 

We study the relationship between spatial distribution 
of counterion and the effective interaction between the 
RLPEs when all counterions are monovalent. For this, we 
consider the same system specifications as in Sec. IIII CI 
A counterion is said to have condensed on a PE if its sep- 
aration from any of the monomers of the PE is less than 
1.25cr. When the PEs come close by, the counterion may 
be shared by both the PEs, but the qualitative nature 
of the angular distribution remains the same. In Fig. [5l 
we show the dependence of P{9) on 9, where P(9)d8 is 
the probability that the condensed counterion is between 
angle 9 and 9 + d9. The data shown are for A = 9.43, for 
which a pronounced minimum in W(d) was observed (see 
Fig. SI). Here, 9 = or 2ir corresponds to a location in 
the same plane and in between the two PEs, while 9 = n 
corresponds to the counterion being located in the same 
plane but away away from the second PE. When the sep- 
aration d between the PEs is much larger than 1.85a, the 
condensed counterions are distributed uniformly around 
the PE. As the PEs approach each other, the distribu- 
tion develops a sharp peak at 9 = 0, and a broad peak 
at 9 = 7r, showing that most condensed counterions are 
located in the same plane as that of the PEs, and mostly 
between the PEs. For separations less than 1.85er, the 
peak at 9 = is shifted to higher values, showing that 
the condensed counterions are expelled out of the plane 
of the two PEs, though they remain in the region be- 
tween the PEs. This result is in accordance with Man- 
ning's theory^ of shared counterions being the origin of 
attraction between two similarly charged PEs. The ex- 
clusion of counterions unscreens the interaction between 
the similarly charged monomers of the two PEs, caus- 
ing the repulsion between them. These observations are 
consistent with the nature of the effective interaction po- 
tential between the PEs, i.e., attraction for intermediate 
separation(1.85er < d < 4er) and repulsion at very small 
separation (d < 1.85a). 

IV. CONCLUSIONS 

We studied, using molecular dynamics simulations, 
the role of the valency of the counterions on the aggre- 
gation of similarly charged RLPE chains. We showed 
that monovalent counterions can mediate the aggrega- 
tion, provided the linear charge density of the PE back- 
bone is larger than a critical value. The absence of ag- 
gregation in the presence of monovalent counterions in 
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FIG. 5. The angular distribution P(6) of condensed monova- 
lent counterions around a given PE for A — 9.43. 



earlier simulations and experiments is probably due to 
the linear charge density being smaller than the criti- 
cal value. The critical linear charge density for aggrega- 
tion decreases with increasing counterion valency, and is 
found to be close to the critical value for the extended- 
collapsed transition of a single flexible PE chain. We 
also find that, for two parallel RLPEs, the angular dis- 
tribution of the condensed counterions changes with the 
separation between the chains. When the effective inter- 
actions are attractive, we find that the condensed counte- 
rions are shared by the RLPEs, similar to the mechanism 
suggested by Manning 18 . 

It would be interesting to study the aggregation tran- 
sition in more detail. From our simulations, it appears 
that the phase transition is first-order in nature. How- 
ever, using molecular dynamics simulations, it is difficult 
to equilibrate the system in the aggregating regime, es- 
pecially if the initial density is low. Hybrid simulations 
using both Monte Carlo and molecular dynamics might 
be useful to study the transitio n 32 ' 35 . 

Another problem of interest is the kinetics of aggrega- 
tion. From our simulations, it appears that the density 
decay is independent of A and valency of the counteri- 
ons. This may be due to fact that the attraction be- 
tween two RLPEs is short-ranged for all values of A and 
valency. When two aggregates come close by, the proba- 
bility that they aggregate may depend on A and valency. 
This, while affecting the prefactor of the power law de- 
cay, does not affect the exponent, suggesting that the 
kinetics is diffusion limited. Understanding the kinetics 
better, and also its dependence on hydrodynamic inter- 
actions will be part of a future study. 

We note that the critical linear charge density for ag- 
gregation of RLPEs in the presence of monovalent coun- 
terions is much higher than that in the presence of mul- 
tivalent counterions. One system which has comparable 
high linear charge densities is solutions of charged worm- 
like micelle o 58 ' 59 . Earlier Monte Carlo simulations of such 
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systems have used high charge densities^. Monovalent 
countcrion induced aggregation, as seen in the current 
paper, may be realized in such charged micellar systems. 
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